#ifndef CSOLVER_H_
#define CSOLVER_H_

#include "Matrix.h"
#include "Vector.h"
#include<cmath>
#include<iostream>

using std::cout;

class CSolver{
public:
	CSolver(){}
	virtual ~CSolver(){}
	Vector solveGauss(Matrix &m,Vector &v){
		int i,j,k;
		int n=m.row();
		Vector r=Vector(n);
		double **a=m.getValues();
		double *b=v.getValues();
		double *x=r.getValues();
		double Sum, p;
		for (k = 1; k <= n - 1; k++)
			for (i = k + 1; i <= n; i++) {
				p = a[i][k] / a[k][k];
				for (j = 1; j <= n; j++)
					a[i][j] -= p * a[k][j];
				b[i] -= p * b[k];
			}
		for (i = n; i >= 1; i--) {
			Sum = 0;
			x[i] = 0;
			for (j = i; j <= n; j++)
				Sum += a[i][j] * x[j];
			x[i] = (b[i] - Sum) / a[i][i];
		}

		return r;
	
	}
	
	Vector solveIterator(const int &max, Matrix &m, Vector &v)
	{
		int i,j,k;
		int n=m.row();
		double error,sum;
		Vector old=Vector::zeros(n);
		Vector r=old;
		for(k=0;k<=max;k++){
			error=0;
			std::cout<<"k= "<<k<<" : ";
			for(i=1;i<=n;i++){
				sum=0.0;
				for(j=1;j<=n;j++)
					if(i!=j)
						sum+=m(i,j)*r(j);
				old(i)=r(i);
				r(i)=(v(i)-sum)/m(i,i);
				std::cout<<r(i)<<"\n";
				error=(error>fabs(r(i)-old(i)))?error:fabs(v(i)-old(i));
			}
			std::cout<<error<<"\n";
		}
		return r;
		
	}

};
#endif
